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Abstract 

This paper presents a model for the simulation of liquid-gas-solid flows by means of the 
lattice Boltzmann method. The approach is built upon previous works for the simulation of 
liquid-solid particle suspensions on the one hand, and on a liquid-gas free surface model on 
the other. We show how the two approaches can be unified by a novel set of dynamic cell 
conversion rules. For evaluation, we concentrate on the rotational stability of non-spherical 
rigid bodies floating on a plane water surface - a classical hydrostatic problem known from 
naval architecture. We show the consistency of our method in this kind of flows and obtain 
convergence towards the ideal solution for the measured heeling stability of a floating box. 



1. Introduction 

Since its establishment the lattice Boltzmann method (LBM) has become a popular 
alternative in the field of complex flow simulations [33] . Its application to particle suspensions 
has been propelled to a significant part by the works of Ladd et al. [22, 23] and Aidun et al. 
[3, 4, 1]. Based on the approach of the so-called momentum exchange method, it is possible to 
calculate the hydromechanical stresses on the surface of fully resolved solid particles directly 
from the lattice Boltzmann boundary treatment. In this paper, the aforementioned fluid- 
solid coupling approach is extended to liquid-gas free surface flows, i.e., the problem of solid 
bodies moving freely within a flow of two immiscible fluids. We use the free surface model 
of [21, 30] to simulate a liquid phase in interaction with a gas by means of a volume of fluid 
approach and a special kinematic free surface boundary condition. I.e., the interface of the 
two phases is assumed sharp enough to be modeled by a locally defined boundary layer. 
This boundary layer is updated dynamically according to the liquid advection by a set of 
cell conversion rules. 

This paper proposes a unification of the update rules of the free surface model with those 
of the particulate flow model, which also requires a dynamical mapping of the respective solid 
boundaries to the lattice Boltzmann grid. As described in [5] , the resulting scheme allows full 
freedom of motion of the solid bodies in the flow, which can be calculated according to rigid 
body physics as in [20]. We demonstrate the consistency of the combined liquid-gas-solid 
method by means of a simple advection test with a floating body in a stratified liquid-gas 
channel flow, and discuss the main source of error in the dynamic boundary handling with 
particles in motion. 
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Figure 1: D3Q19 stencil. The weights are w = 1/3 for C, Wx,..,w e = 1/18 for W, E, N, S, T, B, and 
w 7 , .., w ls = 1/36 for TW, TE, TN, TS, NW, NE, SW, SE, BW, BE, BN, BS. 
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We further apply our method to the problem of rotational stability of rigid floating 
structures. This kind of hydromechanical problems typically emerge in marine engineering, 
where the floating stability of offshore structures is of concern [19, 12], such as the stability of 
a ship in a heeled position. Because of the static nature of the addressed problems, numerical 
issues arising from hydro-dynamic effects can be widely discarded, which makes them well- 
suited for the verification of the involved force calculations. In addition to that they provide 
a possibility to check the convergence of the simulated liquid-gas-solid systems into a state 
of equilibrium. We succeeded in showing that basic convergence is given, provided that 
adequate spatial resolutions are chosen. For the special problem of the floating stability of 
cuboid structures, convergence of numerical simulations towards the analytical model was 
obtained. 

The idea of evaluating the simulated floating stability of rigid bodies was inspired by 
[14], who proposed it as a test case for a Navier-Stokes based simulator originally developed 
for the estimation of "green water" loads on ship decks [13]. As mentioned above, lattice 
Boltzmann based fluid-structure interaction techniques have been developed for the simula- 
tion of particulate flows. We were not able to find publications discussing the application of 
LBM-based free surface flows in interaction with floating structures. To our knowledge, the 
only approach to handle similiar problems is the one proposed by JanBen [18]. 

2. Method 

2.1. Isothermal D3Q19 Lattice BGK Method 

We assume the D3Q19 lattice model for 3-dimensional flows [31], with a set of iV = 19 
discrete lattice velocities q (i = 0, .., N — 1). For the theoretical considerations in this 
section, however, we will often fall back implicitly to the native D2Q9 model, as it simplifies 
explanations and figures if they are 2- dimensional. The lattice velocities c\ (also called lattice 
directions or lattice links) with their respective weights Wi (i — 0, .., N — 1), as shown in Fig 
1, are 

r (0,0,0) 

Ci= I (±1,0,0),(0,±1,0),(0,0,±1) (1) 
[(±1,±1,0),(0,±1,±1),(±1,0,±1). 
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In the following i is used to denote the index of the lattice velocity c- with cj = — q. 

Let s x , s y , s z be positive real numbers divisible by the spatial resolution, constant 8 X . 
The domain [0, s x ] x [0, s y ] x [0, s z ] is divided into cells, i.e., cubic volumes of length S x , which 
yields a computation domain of ^ x |^ x |^ discrete lattice cells. Spatial quantities like s x , 
s y , s z and 5 X are commonly given in a certain unit of length (e.g., metres). However - when 
dealing with LBM specific computations - dimensionless lattice coordinates are used: Spatial 
coordinates are thus given in the following as multiples of 5 X . By speaking of a cell (x, y, z), 
where x, y and z are positive integer numbers we mean the lattice cell with respective volume 
[x, x + 1] x [y, y + 1] x [z,z+ 1] in the lattice. We refer to the point (x + 0.5, y + 0.5, z + 0.5) 
as the cell center. For each lattice direction i = 0, .., N — 1 we name fi(x,t) the particle 
distribution function (PDF) of the direction q in cell x and of time step t. 

The lattice BGK propagation scheme can be derived from the classic Boltzmann equation 
with the collision operator substituted by the Bhatnagar Gross Krook (BGK) operator [27]. 
Including an external force term Fi, the lattice BGK (LBGK) equation reads 

fi(x + 5tcl, t + 6 t ) = fi(x, t) [fi(x, t) + f eq i (p(x, t),u(x, t))] - 5 t Fi. (2) 

T 

t is the dimensionless relaxation time and related to the kinematic viscosity v by 

v + l/2c% 



T 



c 2 A 



The lattice speed of sound c s is a mo del- dependent constant. For the D3Q19 model it is 
c s = The equilibrium function is therewith given as a so-called low Mach number 

expansion of the Maxwell distribution function [17], 



feq,i(p,u) = pWi 



Cj U (Cj u) u u 

1 + ^+ v 



c 2 s 2ct 2cl 
and is valid only for small flow velocities, where the following constraint holds: 



(3) 



u u 



Ma := — < 1. (4) 



c 2 

1 s 

The external force term is used to represent gravitation (expressed as acceleration a) in 
the simulation. It is given by [26] 

fcl-u cl T u\ _ 

The local macroscopic quantities, density p and fluid momentum pu, are obtained as moments 
from the PDFs: 

JV-l 

P = E /» ( 6 ) 

i=0 
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N-l 

pU=Y^ ( 7 ) 

i=0 

N-l 

puu T + pi + S = 22 clcffi (8) 

i=0 

The last equation contains the momentum flux tensor puu T and the stress tensor, where 
S represents the shear stresses of the flow. For a system in equlibrium, i.e. fa = f eq>i 
(i — 0..N — 1), the stress tensor S vanishes. 

Within this lattice Boltzmann model, the pressure is linearly related to the density, by the 
equation 

V = c 2 sP . (9) 

When a simulation includes an external force, e.g. in order to simulate the effect of 
gravity, stability concerns limit the density gradient and therefore the hydrostatic pressure 
gradient. Buick and Greated [7] give the following incompressibility condition which relates 
the external force a to the projected height of the simulation domain in the force direction. 

(l x ,ly,l z )-a<(? 8 . (10) 

l x = j 1 , l y = y- and l z — y are t ne respective numbers of fluid cells cells in the x- , y- and 
z- direction. 

In practice, the calculation according to the LBGK equation (2) is split up into two 
steps. First, a propagation step, the stream step, which propagates the local PDFs of each 
cell along the corresponding lattice link into the neighboring cell. The central one with zero 
lattice velocity remains on the same cell. 

ft(x + c i ,t+l) = f i (x,t) (11) 

Now, in the collide step, the relaxation towards the local equilibrium is performed for each 
cell with the set //, % = 0, .., N — 1 of PDFs. The external force Fi is added. This yields the 
PDFs of the succeeding time step. 

fi(x, t + l) = //(£, t + 1) - i [//(£, t + 1) - f eq ,i(p, u)} - F. (12) 

For a step by step derivation of the lattice Boltzmann method see [16]. 

2.2. Previous Work: Free Surface LBM 

The main characteristic of the free surface model given by [21] is that the dynamics of 
the second phase are neglected. It assumes a system of two immiscible fluids, in which the 
dynamics of the first phase govern the flow completely. We refer to this first phase as the 
liquid phase and to the second phase as the gas phase. The layer where the two phases come 
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Figure 2: 2D-Rcprcscntation of a free liquid-gas boundary by interface cells. The real interface (dashed line) 
is captured by assigning the interface cells their liquid fraction. 

i i i I I I rn ngas 

, □interface 
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in contact - the free surface - is assumed to be thin in relation to the spatial resolution 8 X . 
The separate regions of liquid or gas are on the other hand assumed to be large in comparison 
to 5 X . With this assumption, one can simulate the liquid phase by modeling the free surface 
as a special boundary condition. By means of a volume of fluid approach that introduces 
a dynamic fill level (p(x, t) for each lattice cell, partial filling is allowed for cells at the free 
surface boundary. The fill level is the fraction of cell volume 5^ that is currently filled with 
liquid, such that the liquid mass in a cell may be calculated as 

m(x, t) = (f)(x, t) ■ p(x, t)S x . 

Cells that have a liquid fraction < 1 and > thus form a boundary for the fluid referred 
to as the interface, and need special treatment, which is outlined in the following. We refer 
to those cells as interface cells. Fig. 2 shows the interface layer between liquid and gas for 
a virtual free surface, drawn as a dashed line. During the stream step, the fill level may 
change due to the exchange of PDFs with neighboring cells. The mass balance between two 
neighboring cells, x and x + q, is given by 

if x + el is gas, 
if x + ci is liquid, (13) 
fi(x, t)) if x + cl is interface. 

The interface must further fulfill the characteristic of a free surface boundary. Because 
the gas flow is neglected in the model, only the pressure of the gas is taken into account at 
the interface between the two phases. Apart from that, the interface must freely follow the 
liquid flow, which is provided by a special free surface boundary condition. Assuming that 
an approximated normal vector n(x, t) of the free surface is known for the interface cell at 
position x, the PDFs /$ with c^n < (pointing towards the liquid phase) are set to 

//(£,*+ 1) = f eq Ap G (x),u(x)) + f eq: i(p G (x),u(x)) - f{(x,t + 1) (14) 

during the stream step. Here pc(x) is chosen such that by Eq. 9 it matches the pressure p G 
of the gas phase, which must be given. Indeed, Korner et al. [21] show that this results in a 
zero strain rate tensor for the boundary. This can be seen by substituting the reconstructed 
PDFs from Eq. 14 into the formula for the momentum flux and stress tensor, i.e., the second 





I fl(x + c u t) - fi(x,t) 

\ (0(X, t) + <f)(x + Ci, t)) ■ (fi(x + Ci,t) - 
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order moment (Eq. 8) of the PDFs: 

puu T + pl + S = Cicl T f-(x, t + 1) + 

i,Ci T n>0 

^2 ^ [feq,i(PG, U) + feq,i(pG, ~ f{& t + 1)] 
i, cl T n<0 

= ^ ^ T [feq,i(PG,u) + feq,i(PG,u}] 
i,Ci T n<0 
N-l 

= ^ CiCl T feq,i(pG,u). 
i=0 

Because the last sum resembles a system in equilibrium, it follows S = and p = c\pq = Pg- 
The normal vector n(x) is obtained as approximation from the gradient of fill levels within 
a local neighborhood of the cell. It is also possible to include the effect of surface tension 
directly in the boundary condition above. If a constant surface tension parameter a is given 
and the local curvature of the interface k(x, t) is known, the gas pressure is augmented by 
the Laplace pressure to 

Pg =Pg + 2<tk. 

Pohl [30] proposed a method for 3D-calculations that include surface tension. His ap- 
proach is, to extract the local curvature k from the surface normals in a local neighborhood 
clS 8b second step after the normal computation. 

In order to allow free advection of the free boundary between gas and liquid, the state of 
the interface cells must be tracked carefully. Since the interface cells represent a boundary 
for the LBM scheme, it is necessary to assure a closed layer of interface cells around the 
liquid cells. This is facilitated by a set of cell conversion rules. Thus, no direct state 
transition between gas and liquid is allowed. If an interface cell reaches a fill level 4>{x) < —e 
or <p{x) > 1 + e it is converted into a gas cell or a liquid cell, respectively. Since such a 
conversion could introduce holes in the interface, further cell conversions from either gas or 
liquid into interface cells may be triggered. The newly generated interface cells are initialized 
with a fill level <fi = or <p = 1, respectively. This cell conversion scheme has to be extended 
for the incorporation of floating objects, as shown in Sec. 2.4. 

2.3. Previous Work: Particulate Flow LBM 

Walls and other solid obstacles that block the fluid motion must be handled with ap- 
propriate boundary conditions. For this paper we assume that all solid boundaries can be 
modeled by a no-slip condition, i.e. the relative motion of the liquid at the boundary is zero. 
This can be achieved by reflecting distribution functions, that are hitting an obstacle cell 
during the stream step, to the opposed direction. If fi(x) is given with an obstacle cell at 
x + cj it is bounced back according to 

2 

f-(x, t + 1) = fi(x, t) + —WiC^u w p, (15) 
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Figure 3: A box-shaped rigid body mapped into the lattice. With a free surface flow, four different cell types 
have to be distinguished: Liquid, interface, gas and obstacle. 

□ gas 

□ interface 

□ liquid 
^obstacle 




where u w is the velocity of the obstacle. The second term on the right hand side accelerates 
the fluid according to the movement of the wall [22]. If the solid is not moving (u w = 0), 
this term will vanish, and the boundary handling is reduced to a pure reflection. 

From the boundary condition of Eq. 15 one can directly obtain the stresses exerted on the 
boundary. This is known as momentum exchange [22, 23, 28] in literature. The momentum 
transferred to the wall must equal the change of momentum which results from the reflection 
of the PDF, since elastic collisions are assumed. Thus, for a fj reflected to // by Eq. 15, the 
change in momentum is 

Aji = f-cl - fi&i = 2ci(fi + \wicl T u w ). (16) 

c 2 s 

The simulation of immersed rigid bodies behaving according to Newtonian physics is facil- 
itated by discretizing the shape of a body to the lattice. For every time step the no-slip 
boundary is set according to the object position and orientation. If the center of a cell is a 
point within the volume represented by the shape of the body, it is treated as an obstacle cell 
(Fig. 3). The movement of the rigid bodies resulting from the stresses, Eq. 16, including the 
resolution of body-body collisions, is realized in practice by coupling the lattice Boltzmann 
algorithm to a rigid body physics engine as in [20, 11]. The resulting body motion is then 
fed back to the fluid simulation by setting the boundary condition according to Eq. 15 in 
the cells covered by the obstacle. Hereby u" w is given by the velocity of the body. 

A given particle must behave according to the stresses exerted by the fluid. The net 
force and torque on the particle are obtained from the momentum flux according to Eq. 16, 
summed over the discretized particle surface. Let N be the set of grid points x next to the 
discretized particle surface. For a cell x G N a , let I {x) be the set of indices i G I (x), where 
Aji(x) is defined due to a reflection of PDFs as stated above. Then the discrete surface 
integrals for force F net and torque T, respectively, are given by 

Fnet = E E A ^( f ) • J> ( 17 ) 

x€N a i£l (x) t 

f = E E (z-olxA/^A (18) 

x&N iel (x) 1 
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Figure 4: Allowed state transitions for cell types. With moving obstacles the number of transitions increases 
from four to ten. The free surface type conversions are framed by a dotted line. The dashed arrows indicate 
critical state changes from solid to fluid. 



In the latter sum, o is the center of gravity of the object. A non-zero torque acting on a 
freely moving rigid body will induce rotational movement of the body. 

The momentum exchange method has become a kind of quasi-standard for the simulation 
of particle suspensions with Lattice Boltzmann and it has been mentioned in a large number 
of publications with only slight modifications. The main difference here to the works of Ladd 
[22, 23] or Aidun et al. [3, 4] is, that in our case the obstacle cells are always fluid blocking 
nodes. In the mentioned references there are "virtual" fluid nodes inside the discretized 
object shapes, which can be reactivated when uncovered by from the object. This approach, 
however, would be difficult to maintain in the presence of a liquid-gas flow, since it is unclear 
how the liquid-gas interface would have to be treated when intersecting with the volume of 
a particle. 

24. A Liquid- Gas- Solid LBM 

With the free surface model of Sec. 2.2 and the additional presence of moving obstacles 
from Sec 2.3, there are now two different boundaries for the liquid phase that must be treated 
dynamically. Fig. 3 illustrates the situation for a box floating on the free surface of a liquid, 
simulated according to Sec. 2.2. Since there are now four different types of cell states to 
consider (liquid, interface, gas, and solid) that may change according to fluid flow and object 
movement, the need for a sophisticated cell conversion algorithm arises. In the following, 
we show that these conversions can be organized in such a way that a consistent boundary 
around the liquid phase is assured in every time step. 

Fig. 4 shows that the total number of conversions consists of ten state transitions. The 
transitions at the bottom are those of the advection scheme of the free surface model of Sec. 
2.2 and need no further examination. The remaining transitions arise from the mapping of 
the solid to the grid. We assume the velocity of the object is small enough for conversions 
into obstacle to occur only at positions in the neighborhood of other obstacle cells; while 
conversions from obstacle to one of the three fluid states (liquid, interface or gas) can occur 
only at positions in the neighborhood of fluid cells. This is no additional limitation, but 
a direct consequence of the low Mach number limitation of the lattice Boltzmann scheme. 
Consider now an object penetrating the free surface as shown in Fig. 3. While conversions 
into obstacle cells are generally safe, the conversions from obstacle into fluid could lead to 
invalid lattice configurations with holes in the interface layer. Thus the correct cell type, 
i.e., liquid, interface or gas, needs to be determined such that a valid lattice configuration 
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is adhered for the next time step. This is achieved by an update rule based on the local 
neighborhood of a given obstacle cell. 

We define the non-obstacle neighborhood of a cell as B(x) := {x+q | x+q is no obstacle}. 
For a prior obstacle cell x to be converted into fluid, the state is determined as follows [5]. 

• B(x) contains no liquid: x is converted into gas. 

• B(x) contains no gas: x is converted to liquid, p(x) is interpolated from B(x), and 
the fi(x) are set to f eq ,i(p(x) , u w ) , where u w is the velocity of the point on the object's 
surface, which is closest to the cell center of x. 

• B{x) contains liquid and gas: x is converted into interface. p{x) is interpolated from 
the non-gas cells in B(x), and the fi(x) are set as in the preceding case. Additionally 
a fill level 4>(x) is chosen by interpolation of the fill levels of the interface cells in B(x). 
It is also possible to include the gas cells with the pressure pc in the interpolation of 
p, since pressure can be seen as a continuous quantity across the free surface boundary 
ofEq. 14. 

With the free surface model of Sec. 2.2, the gas phase is only taken into account in terms 
of a pressure force exerted onto the free surface. This pressure forces, however, are also 
acting on the surface of the solid bodies. So, in analogy to the boundary condition for the 
free surface, Eq. 14, we use the equilibrium function, Eq. 3, to construct PDFs f eq ,i(pG, 0) 
for the lattice links from gas, or interface cells to the obstacle cell of an object. From these 
the momentum transferred to the body can be calculated from Eq. 16. The value of the 
constructed PDF is linearly interpolated with the liquid's distribution function Thus, in 
presence of gas in the neighborhood of an obstacle cell, the momentum acting on the body 
locally is given by 

Al = 2&i U ■ h + (1 - <f>) ■ feqjipa, 0) + ^ Wi C^U^j . (19) 

For a gas cell, = 0, only the equilibrium function of Eq. 16 remains, all other parts vanish. 

3. Simulation of Floating Bodies 

In the following sections we present the results of calculations performed with the method 
of Sec. 2.4, i.e., various simulations of a rigid body exposed to a free surface flow. Since 
many details of the method are owed to previous particulate flow simulations, we briefly 
recapitulate the most relevant results in Sec. 3.1. We then focus on the evaluation of 
forces and motion of partially immersed rigid bodies, Sec. 3.2, in order to prove the basic 
consistency of the method. When dealing with floating objects in gravity-stratified free 
surface flows, the buoyancy forces on the discretized objects turned out to be one of the 
main sources of error. In the last part, Sec. 3.3, we apply our method to the problem of 
floating stability of rigid bodies. 
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3.1. Force Computations with Spherical Particles and Staircase Approximation 

It has to be emphasized that due to the discretization of the body shapes to lattice sites as 
described in the previous section, a certain a priori error is introduced. This is a well-known 
issue in LBM-based particulate flows, and there has been a lot of research regarding the 
treatment of complex, curved boundaries in lattice Boltzmann over the last decades [31, 2]. 
These schemes allow second order spatial accuracy in the treatment of curved boundaries. 
Most of them are directly applicable for moving boundaries, such as particles. However, in 
the kind of flows which we are going to discuss in the final part of the paper, interpolation 
based schemes like [6], or [35], which have been applied with success to particle simulations 
in the past, turned out to be little useful. 

The solid particle model described in Sec. 2.3 which we use throughout this paper is 
identical to the one used for drag force computations on suspended particle agglomerates 
in [11, 20], where the net force deviation was studied with spherical particles in stokes 
flow. The calculations were performed with various degrees of spatial resolution and a set 
of different lattice relaxation parameters r in [11] with a sphere fixed to the center of a 
channel flow. The results were found to be in good agreement with theoretical expectations. 
Second order boundary conditions surpassed the staircase approximation under all tested 
circumstances. In [20], sinking sphere simulations were conducted in order to evaluate the 
method with moving particles as proposed in [1]. These simulations included the influence of 
gravity, and used the unified boundary treatment by [35] for the solid boundaries: A particle 
of density p s > 1 accelerates until the frictional force imposed by the fluid surrounding 
balances its weight and a terminal sinking velocity is reached. The results showed small 
periodic oscillations in the measured force on a sphere that moves with a quasi- const ant 
velocity through the lattice. This error can be attributed to the staircase approximation of 
the boundary of the body. Fig. 5 shows the volume error for a sphere of radius 5 5 X when 
traveling along an axis-aligned path. The number of cells actually covered by the sphere 
is depending on its exact position relative to the lattice. With the discretized shape of a 
particle changing, the set of lattice links included in the discrete surface integral around its 
boundary changes and so does the net force, Eq. 17, and torque, Eq. 18, respectively. The 
interpolation based schemes applied in [11, 20] however, approximate a curved boundary by 
calculating its exact distance from the cell center along the lattice link of reflection. This 
information is used to correct the reflected PDFs from Eq. 15. Hence, the error arising in 
the buoyancy force from changes in the discrete volume of an object cannot be corrected 
by purely local refinements. For liqud-solid suspensions without a free surface the problem 
of buoyancy force fluctuations with fluid-blocking obstacles is usually negligible, or can be 
corrected by workarounds as reported in [15]. The following subsections discuss this issue 
for the free surface case. 

3.2. Forces on Buoyant Rigid Bodies 

Due to the staircase approximation of rigid body shapes, not only is the buoyancy force 
depending on the object's position, but - for partially immersed objects - a staircase function 
with respect to draft. This can be best explained for an axis-aligned cube of side length s 
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Figure 5: Discretized volume for a spherical particle traveling along an axis- aligned path with a velocity of 
1CP 4 . Ideal volume: 523.6. The graph repeats periodically every 10000 time steps 



in a liquid at rest with a constant water level. The buoyant force is then given as a linear 
function with respect to the draft d: 

Fb = pgs 2 d. 

Because of the discretization error, d is only exact up to the lattice resolution 5 X . This can 
keep buoyant objects from coming to rest at their equilibrium position. Especially for objects 
with lattice-aligned planar faces such as a cube this would sometimes cause an unphysical 
"wobbling" behavior with the object jumping between two points slightly below and above 
the ideal floating position. If on the contrary a fixed cuboid is considered and the gauge of 
the liquid is varied immersing more or less of the object's surface, the changes in the buoyant 
force on the object would be continuous. This is because the volume of fluid-based interface 
treatment is not limited to the spatial accuracy of the lattice. The fill level is taken into 
account in the force calculations in the interface cells at the triple line of liquid, gas and 
solid by Eq. 19. Test cases with a linearly rising water level around a partially immersed, 
axis-aligned, fixed cuboid led to a linear increase in the buoyant force. 

3.2.1. Free Advection Test 

To test the consistency of the unified liquid-gas-solid algorithm, a free advection test was 
performed in a channel flow with a constant homogeneous velocity in the laminar regime 
along the x-direction, and a gravitational field along the z-axis. A spherical particle of 
density p s = 0.5 was placed freely floating in the center of the channel to be accelerated by 
the hydrodynamic forces acting on its surface. Fig. 6 gives an impression of the simulated 
setup. Because the gas phase is neglected in the model it must not hinder the particle motion 
in any way. This means, that in both cases the final velocity of the particle must equal the 
fluid velocity. After the systematic drag force evaluation cited in the preceeding section, 
this test might seem trivial at the first spot. But in the presence of a density gradient and 
because of the changes in the discretized particle shape (Sec. 3.1), there are at least two 
important sources of error to consider. Hence, in the following we also examine the results 
at critical flow velocities that have the same order of magnitude as the expected errors. 

The channel size measured 60 x 40 x 30 lattice units, and the particle diameter was 
12S X . All long sides of the channel were no-slip moving walls (Eq. 15), in order to obtain 
a homogeneous velocity profile across the domain with a constant flow velocity u along the 



11 



Figure 6: Schematic view for a free advection test with a floating sphere in laminar free surface flow. The 
channel is driven by a moving wall boundary condition at the bottom and its longitudinal side planes. The 
channel entry and outflow are connected with each other by means of a periodic boundary condition. 
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x-axis. The left (x — 0) - and right (x = s x ) - boundaries were connected periodically, by 
copying PDFs streamed outside at one end into the cells at the opposing end of the domain. 

The whole process of testing was performed at four different values of T1..4, namely 0.62, 
0.8, 1.1, 1.7 a gravitational constant of g = 10 -5 , and four different channel velocities M1..4 = 
10~ 4 , 5 • 10~ 5 , 2.5 • 10~ 5 , 1 • 10~ 5 . In order to check for possible sources of errors the scenario 
was run in a sequence of increasing complexity, starting with a simplified version of the 
simulation experiment. The first series were the following. 

• Channel fully filled with liquid without gas; without gravitational field. 

• Repetition without gas, but with gravitational field as stated above. The particle 
density was set to 1.0 in this case. 

In the first run the particle velocity matched the liquid velocity to the last digit without any 
deviations. For the second run we were expecting errors due to oscillations in the buoyant 
force as described above. These were clearly measurable, and largest for the lowest lattice 
viscosity T\ = 0.62, where the z-component of the object velocity oscillated within the interval 
[—g, +g\. For larger r values the effect would be damped (e.g., by a factor of 0.1 in the case 
of T4). There were no significant deviations in the other velocity components of the particle. 

The next test series consisted of the channel partially filled with liquid with a planar free 
surface at z = 19.5 under the influence of gravity. The particle was removed to check the 
behavior of the interface cells. The free surface remained planar under all tested circum- 
stances, as expected. For the lowest value of r = 0.62, locally spurious interface velocities 
were occuring. However, the absolute deviation from the bulk viscosity was always below 
10~ 5 , restricted locally to parts of the interfacial layer, and therewith only critical for low 
velocities. 

We finally tested the advection of a solid body in the free surface flow by placing the 
particle with density p s = 0.5 inside the channel. The fluctuations in the lifting force were 
comparable to those measured in the all-liquid domain under gravity with all tested cases. 
However, here the error also became visible in the in the other components of the net force 
on the object, and thus introduced an error in the streaming velocity of the particle. We 
believe that this is caused by the errors introduced by the refilling step of the algorithm, 
when obstacle-to-fluid conversions occur for an object sticking through the interface layer. 
As described in Sec. 2.4, we initialize these cells to a local equilibrium with parameters 
interpolated from a local neighborhood. This is likely to introduce a certain error that cannot 
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be balanced in the asymmetrical case when an obstacle is at the boundary between gas and 
liquid. Even in case of liquid-solid flows without a gas phase, the problem of initializing new 
liquid cells is non-trivial [8]. 

We have not evaluated those cases, where the empty channel showed spurious currents. 
These were: T\ = 0.62 with velocities 1*2, •■,1*3. In the remaining cases the sphere would 
follow the current with the velocities given in Tab. 1. For higher flow velocities (or lower 
values of g), this error would no longer be critical, if even apparent in the sphere velocities. 

Tabic 1: Errors in the velocity of the sphere for lattice gravity of g = 10~ 5 : Critical if r is low and the 
characteristic flow velocity of the same order as g. 



T 


Ml 


u 2 


u 3 


W4 


0.62 


l.le-4 


n.e. 


n.e. 


n.e. 


0.8 


1.05e-4 


5.5e-5 


3.0e-5 


2.0e-5 


1.1 


1.01e-4 


5.15e-5 


2.8e-5 


1.3e-5 


1.7 


l.e-4 


5.05e-5 


2.55e-5 


l.le-5 



The the systematic error in the staircase approximated objects is negligible as long as the 
characteristic velocity U of the flow is large in comparison with the gravitational constant. 
However, it should be considered in the Froude number scaling [19] (Fr = U/y/gL, for a 
characteristic length L) of simulations. On the other hand the question arises, whether a 
more accurate treatment of particle boundaries can be achieved that does not suffer from 
the presence of a hydrostatic density gradient. We were unable to find any literature that 
reports the problem explicitly. Besides increasing the spatial resolution, or applying local 
grid refinement methods that entail complex restructuring of the whole implementation of 
the method, we found at least one approach that could provide a significant improvement. 
Based originally on Chen et al., a volumetric interpretation [9] of the lattice Boltzmann 
method allows the direct calculation of the momentum flux onto finite surface elements. 
This can be exploited to realize boundary conditions, like [10, 34, 32], that offer sub-grid 
scale accuracy. [32] presented a variant that allows an exact description of the surface of 
obstacles including surface stress calculation which is independent of surface motion. As 
the free surface boundary treatment, Eq 14, is in principle also a volumetric approach [21], 
this would yield a further unification of the method described in this paper. We expect to 
explore this issue in future publications. 

3. 3. Stability of Floating Bodies 

We now apply the method on the problem of floating stability of non-spherical rigid 
bodies. Archimedean law appears to be quite intuitive, which states that the buoyancy 
experienced by an immersed object is equal to the weight of the liquid displaced by V, the 
immersed volume part of the body (i.e., the part below the water plane). According to 
that rule, a body with a specific material density p s < 1.0 will sink down or rise up until 
buoyant force Fb = pgV and gravitational force Fq balance each other vertically. In practice, 
however, even if this balance of forces is given, the resulting equilibrium is often found to be 
unstable. Consider for instance a cube in 2-D as shown in Fig. 7. In both positions exactly 
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Figure 7: For both positions of the cube of specific density p s = 1/2 hydrostatic equilibrium is given. 
However the equilibrium on the left hand side is unstable. A stable equilibrium can be found with the cube 
rotated 45 degrees. 
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Figure 8: Stable floating positions for a cube of material density p s = 1/4 (left) and p s = 3/4 (right), 
respectively. Each time one corner of the cube is in place with the water line while another one is exactly 
half-immersed. 




Figure 9: Unstable equilibrium of a cube. Because of the slight horizontal displacement of the center of 
buoyancy B relative to the center of gravity of the cube G, a moment occurs that pushes the cube further 
from its upright position. The center of buoyancy B is constructed as the center of gravity of the trapezoid, 
which the immersed part of the cube makes up together with the fluid surface. To construct B, the base line 
segment u is elongated in one direction by the parallel opposite segment v, while the latter is again elongated 
in the opposite direction by base line u. Next, the shifted end points of each line segment are connected. 
The intersection of the connecting line with the line through the midpoints of u and v gives its center of 
gravity of the trapezoid, B. 
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Figure 10: Box rotating about its longitudinal axis (rolling motion). 
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half of the body volume is immersed and thus the body floats in hydrostatic equilibrium. 
However, it can be shown that only the right hand side with an heel angle of 45 degrees is 
stable. To see this, consider the cube slightly rotated from its upright position, such that 
its central axis (dashed line in Fig. 7) is heeled by a small angle a. The center of gravity G 
of the cube is the point on which the gravitational force Fq acts vertically downward. The 
center of buoyancy B, which is defined as the center of gravity of the immersed part V of the 
body, is slightly displaced to the left as shown in Fig 9. Thus with Fb acting on B vertically 
upward, a rotational moment in the heeling direction occurs, pushing the cube further from 
its upright position. In the same way it can be shown that, starting from the 45 degree 
rotated position (right in Fig 7), small angles of heel from this position lead to a righting 
moment in the opposite direction. The result is a stable equilibrium at this position. Fig. 8 
shows two more stable floating positions for the cube densities p s = 1/4 and p s = 3/4. 

3.3.1. Equilibrium Position of Floating Cubes 

Following [14], we check the rotational equilibrium of floating cubes of density p s = 0.25, 
0.5, and 0.75 in simulations. With the given densities 0.25 and 0.75 a stable angle of heel is 
found at 26.565° from the upright axis of the cube as shown in Fig. 8. The cube of density 
0.5 should calibrate to 45°. Since we are dealing with 3-dimensional objects, the 2D cube 
of side length s is represented by a cuboid of dimensions 2s x s x s, that will rotate around 
its longitudinal axis only, as in the schematic of Fig. 10. Inside a domain of 130 x 40 x 40 
cells, the lower half was filled liquid with lattice viscosity r = 1/1.9 = 0.526 and initialized 
with hydrostatic pressure (g = 7.5 ■ 10~ 4 ). Three boxes with s = 16 were placed inside 
the basin, slightly heeled from their axis-aligned upright position by an angle of 2.86°. The 
simulation was run until all fluid motions ceased and the objects came to rest. Fig. 11 shows 
visualizations of the simulated scenario at various time steps. Since the initial positions of the 
objects were far from the equilibrium state, their motions would cause a wavy disturbance 
of the liquid in the basin. It can be clearly seen from these visualizations that even at this 
low resolution a strong convergence towards the ideal equilibrium was given. 

3.3.2. Righting Stability Moment of Wall-Sided Structures 

Consider again a floating structure heeled about a certain angle a about its upright axis. 
The intersection of the line with direction Fb through B with the upright axis of a body 
is defined as the metacenter M. If M lies above the center of gravity G then the floating 
position of the body will be stable, otherwise unstable (see Fig. 9). This concept originates 
from marine engineering and is commonly used in the characterization the floating stability 
of ships and other offshore structures. A comprehensive discussion can be found in [19]. The 
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Figure 11: Box objects of density 0.75, 0.5 and 0.25, respectively, striving towards vertical and rotational 
equilibrium. The first three pictures show the system in motion. The last picture was taken after all visible 
motions had ceased. All snapshots were visualized from simulation output with the ray tracer [29]. 
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lever arm of the torque arising from the horizontal displacement of the mass centers B and G 
is given by GMsina, and hence the righting stability moment of the body can be calculated 

as 

m s = Fb ■ GM sin a. 

The curve m s (a) = Fb ■ CM sin a is called the stability curve of the floating object. This 
curve can be derived analytically for many relevant situations, depending on the shape of the 
immersed structure. An important case is that of a wall-sided structure, where the points 
of the structure's surface sinking under or rising above the water plane upon the considered 
angles of heel are parallel opposing sides when the structure is upright. If a wall-sided 
structure is given, then, according to [19, 12], the Scribanti formula states 

B M = — ■(! + - tan 2 a), 

V 2 ; 

where Bq is the center of buoyancy of the upright position, I is the moment of inertia of the 
water plane and V stands for the immersed volume of the whole body. In case of a cuboid 
of length /, width b, height h and draft d (— p s • h - obtained from the non-heeled position), 
the formula can be simplified by putting 

I_ 

V ~ 127' 

and with K as the keel point of the cuboid, an expression for GM can be introduced as 

gm = k~w q + B^M - KG. 

From this, the ideal stability curve of an arbitrary cuboid can be calculated. Fig. 12 
shows the schematic stability curves for some cuboid structures of density p s = 1/2 and 
various width:height ratio. In case of a cube, the stability curve has a negative slope at the 
origin. This shows, that the upright position is an unstable equilibrium. Increased width 
yields a higher floating stability. 

To further test the force calculation on floating bodies, the righting stability moment of 
a floating cube at different angles of heel was measured in lattice Boltzmann simulations 
and compared to the ideal stability curve of the structure. In all of the following cases, the 
gravitational constant was chosen to be g = 10 -4 and a partially filled basin was initialized 
with a pressure gradient according to hydrostatic equilibrium. The relaxation time was 
t = 1.0. The initial width to height ratio of the box was 6 : 4 (a), and then in a second pass 
a lower ratio of 5 : 4 (b) was chosen. Each time, the error was examined at three different 
resolutions, i.e., the box size (b x h) in lattice units was 24 x 16, 48 x 32, 96 x 64 and 20 x 16, 
40 x 32, 80 x 64, respectively. The length side along the axis of rotation was always chosen as 
I = 2h. The box was positioned axis-aligned with the lattice, exactly half-immersed relative 
to the free surface plane. During the simulation, the box was fixed to a constant position 
and the angle of heel of the half-immersed box was varied with a in 0°..30° around the 
longitudinal axis of the box through its center of gravity. 
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Figure 12: Stability curves for various floating cuboids of density 1/2. From the negative slope of the graph 
for the cube case (width:hcight ratio 4 : 4) at 0° angle of heel, it can be deduced that the upright position 
is unstable. The stability increases if the width is enlarged. 




Figure 13: Stability curve of a heeled cuboid at different resolutions dx = 1, 0.5, 0.25. The legend reads the 
width of the box in lattice units at the different resolutions. There are significantly less deviations in the 
torque of the more stable box of width per height ratio 6:3 (a) compared to the box shape of ratio 5:4 (b). 
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Fig. 13 shows the simulation results compared to the ideal stability curve. The error was 
much larger in the second case (b), since the stability lever arm of the box is shorter at the 
lower width per height ratio. Thus even small errors in the box discretization have a visible 
influence on the behavior of the box. However, in all cases it can be seen that the errors are 
decreasing at higher resolutions. Convergence towards the ideal values was clearly given. 

4. Conclusion 

In this paper we described a method for the simulation of liquid-gas-solid flows by means 
of the lattice Boltzmann method. Based on previous works, the momentum exchange method 
is used to obtain fluid stresses at the surface of rigid bodies, which allows to calculate 
the resulting net force and torque. A novel set of cell conversion rules (Sec. 2.4) allows 
unrestricted movement of the solid bodies within the domain, including penetrations of the 
free surface. We demonstrated the consistency of the method in a free advection test with a 
particle following a homogeneous free surface channel flow under the influence of gravity in 
Sec. 3.2. 

For further validation of force and torque calculations on non-spherical rigid bodies, 
we applied the method to the classic mechanical problem of floating stability of immersed 
structures. Basic convergence towards the equlibrium state was successfully checked for 
the case of box objects with square sides. We took advantage of the fundamental stability 
formulae for wall-sided structures and validated the righting moment of box shaped objects 
under varying angles heel. The results were found to converge towards the ideal values with 
increased spatial resolution. 

As a next step it would be interesting to move on from these hydrostatic stability examples 
towards dynamic ones, that include coupled interaction of liquid and solid objects. We will 
present further results together with method improvements in future publications. 
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